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Abstract 

Using the recently proposed density matrix renormalization group tech- 
nique |l| we show that the magnons in the S = 1 antiferromagnetic Heisenberg 
chain effectively behaves as bosons that condense at a critical field he- We de- 
termine the spin-wave velocity, v = 2.49(1), as well as the gap A = 0.4107(1) J. 
75.10.-b, 75.10.Jm, 75.40.Mg 
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It is by now well established both experimentally p[ and theoretically 1^-0 that the 
5 = 1 antiferromagnetic chain has a gap, A, to a triplet excitation above a singlet ground- 
state. Thus the magnetization, M, remains strictly zero up to a critical field he = A. For 
the quasi-one-dimensional system NENP the critical field is found to be about lOT for 
a field applied along the symmetry axis. The behavior of M{h) just above he has been the 



subject of some theoretical work |9HTl|. The existing experimental results appear to 
be dominated by extraneous effects such as an off-diagonal alternating component of the 
gyromagnetic tensor [|13|, inter-chain coupling, anisotropies and impurities. 



This problem can be solved using a boson quasi-particle model involving triplet magnons 



with repulsive interactions for parallel spins |TT|. This model predicts that h^ = A. (We set 



the Bohr magneton and ^f-factors to one.) At this field the "rest-mass energy" of the magnons 
is exactly cancelled by their Zeeman energy and one- dimensional Bose condensation occurs. 
What prevents a catastrophe from occurring at he is the inter-magnon repulsion which 
governs the behavior of M{h) above h^. To calculate M{h) we need to calculate the energy, 
E{M) of a very dilute system of M polarized magnons in their ground-state. M{h) is then 
found by inverting the equation: h = dE/dM. It was argued in Ref. that the behavior 
of E{M) is the same as for a system of non-relativistic non-interacting fermions: 

E=iA-h)M + Lr^'-^. (1) 

Here v is the velocity, determined from the single magnon dispersion relation at low energies: 

uj{k) = A + v\k - 7r)V2A + 0{[k - 7r]3), (2) 

L is the length of the system and kp is determined from the particle number: 

/I^F dk 
— = Lkp/Ti. (3) 
-kp 271 

This gives, E = M(A - h) + {v-nfM^/QAL'^ + O(M^) and hence 



M/L = ^{h - A)2A/7rf , (4) 
up to terms of higher order in h — A. 



This formula was first given based on a non-interacting fermion model of magnons 0. It 
was later argued to also arise from an interacting bosonmodel sufficiently close to h^, and 
hence to be exact. It is expected to be valid for very general short-range repulsive interactions 
between the spin-polarized magnons. It should hold for arbitrarily weak interactions for 
long enough chains. The reason is that when the average inter-magnon spacing is very 
large compared to the range of the repulsive interaction the multi-magnon wave-function, 
'^m{xi,X2, ■■■Xm), can be approximated by a free fermion (Bloch) wave-function multiplied 
by the sign function e(xi,X2, ■■■Xm) which has the value ±1 and changes sign whenever two 
particles are interchanged of the following form 

^M(a:i,a;2,...a;M) = ^=e(xi,X2, ...a;Af) Yj Y[^iixiJ^2ixi^) ■ ■ ■^M{xij„)sgio.P. (5) 

Here P denotes the permutation and sgnP the sign of the permutation. ipi{x) is a single- 
particle non-interacting wave-function depending on the wave- vector ki. is symmetric 
as required by Bose statistics, is a solution of the non-interacting Schroedinger equation 
almost everywhere (ie. except where two or more particles intersect) and vanishes whenever 
two or more particles come together. As such it is expected to become exact in the dilute 
limit and hence to give exactly the magnetization as /i — /ic. 

Until very recently it has been essentially impossible to test this prediction by numerical 
simulations. The reason is that extremely long systems are required. To simulate the M- 
magnon problem in the dilute regime we need a length L » where ^ 6 is the 
correlation length or the approximate range of the inter-magnon interaction. In practice 
it appears that L > 30M is required. Before last year, the longest chains that had been 
studied accurately had L = 32 so good results were only available for a single magnon and 
it was impossible to study magnon interactions. This situation has now completely changed 
thanks to a breakthrough in the real-space renormalization group technique made by Steven 
White which makes it quite feasible to study chains of length 100 or longer using a density 
matrix renormalization group (DMRG) approach. We will present results here on chains 
of length up to 100 containing up to three magnons. White and Huse obtained related 
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numerical results independently. We analyze our results in a different way which establishes 
Eq. (|^). Our results indicate quite convincingly that the lowest energy two or three magnon 
state has the form discussed above; namely the lowest-energy free fermion wave-function 
multiplied by the sign function. We establish this result in two ways. Firstly we study 
S^{x) = Y^jLi — Xj), where Xj is the position operator for the j'th particle, showing that 
it has the expected form: 

< S"'(x) >= M j dx2dx'i...dxM\'^ m{,x,X2jXzj ■..xm)\^ ■ (6) 

Secondly, we study the finite-size dependence of the energy of the multi-magnon ground- 
state, showing that it behaves as: 

M 

E{M)^Y.^{h). (7) 

i=l 

Holding M fixed and taking L oo, the fcj's are 0(1/L) so the above formula gives 
E(M,L) = M(A — h) + a{M)/L'^. The corrections to Eq. from inter-magnon inter- 
actions are expected to be 0{L^^). This power of L can be most easily understood in the 
M = 2 case. The wave-function contains a normalization factor but is 0{L^'^) over 
the entire region where the interaction is non-negligible since it vanishes proportional to 
k{xi — X2) with k of 0{L~~^). The power arises from squaring the wave-function and 
picking up a factor of L from the integral over the center-of-mass position. We verify that 
the corrections to Eq. are indeed of this order. 

The DMRG method does create certain technical problems because of the fact that 
implementing it optimally requires studying chains with open, rather than periodic boundary 
conditions. Such chains have S = 1/2 excitations localized near the ends which have been 
the subject of a number of experimental and theoretical studies |T^. For our purposes 
they are just a minor annoyance. They are also found to have repulsive interactions with 
the bulk magnons for parallel spins. Thus the magnon wave-functions essentially obey 
vanishing boundary conditions at the ends, with additional corrections of 0{L^^) to the 
energy, Eq. (|^), from the magnon-end excitation interactions. We also considered more 
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general boundary conditions for the magnon wave-functions. These were found to change 
the energy only to 0{L~^) since they lead to negligible changes in the wave-function. 

The DMRG method for open chains leaves only two good quantum numbers the total 
component, S^, and the parity, P. These are conserved under iteration and it is therefore 
possible to work within a subspace defined by these two quantum numbers. We need to 
determine the parity for low-lying states with a given S^. We shall only be concerned 
with chains of even length. For these chains the ground-state is a singlet with even parity, 
0^. Above the ground-state is an exponentially low-lying triplet, 1~. In the thermodynamic 
limit the triplet and the singlet become degenerate and the ground-state four-fold degenerate. 
This spectrum can be seen to arise from the two S = 1/2 end-excitations forming either an 
odd parity singlet or an even parity triplet, in addition to an overall parity flip coming from 
the rest of the ground-state. This parity-flip can be understood from the valence bond solid 
state where we draw two valence bonds emanating from each site. These valence bonds 
represent singlet contractions of pairs of 5 = 1/2's so they have a directionality associated 
with them. When we make a parity transformation we flip the orientation of an odd number 
of valence bonds resulting in a (— ) sign. Thus, the parity, Pe, of a state with no magnons 
present is (+) if the end excitations combine into the singlet and (— ) for the triplet. The 
parity of higher excited states, containing one or more magnons, is a product of three factors, 
PEPswPm- Pm coutalus a contribution of (— ) from each magnon present. This is because 
the magnons are created and annihilated by the staggered magnetization operator, and this 
changes sign upon switching even and odd sublattices. Psw is the parity of the spatial 
wave-functions of the magnons. For instance, for a single magnon, the wave-functions, ipi, 
in Eq. (^) are ipi = ^J2/{L — 1) sin kiX, ki = nrii/^L — 1), with rii odd for even parity and rii 
even for odd parity. We take < x < L — 1 and parity will therefore take x into L — 1 — x. 

For a chain with open boundary conditions the lowest lying state of a given magnetiza- 
tion, M, will have M = m + l, where m is the number of magnons present and the additional 
term, 1, corresponds to the end excitations forming a triplet. In order to minimize the inter- 
magnon repulsion, the wave-function for large L takes the Bloch form of Eq. (J^), with ipi as 
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above and rii = i, in order to satisfy vanishing boundary conditions. Thus Eq. (|l]) becomes 
E„,^,iL) - E,{L) = (A - h)m + ^^^^ ^ g + 0{L-^). = ^. (8) 

We test this formula below for m = 1, 2 and 3. 

The 1~ state becomes degenerate with the ground-state in the thermodynamic limit 
and we shall therefore view this state as the reference state and calculate energy gaps with 
respect to this state and not O"*", as already implied in Eq. (||). We have calculated the gap as 
a function of chain length between this state and three of the low lying states using density 
matrices of the size 243 x 243 keeping 81 eigenvectors of these matrices at each iterations. 
For each of these states we have also calculated < > and < Sj ■ Sj+i > along a L = 100 
site chain using a finite lattice method |jl|. For a discussion of the numerical procedure we 
refer the reader to Ref . ||l[ . A brief summary of our results is shown in Table |. 

The lowest lying M = m + 1 = 2 state corresponds to a state with the end excita- 
tions in the 1~ state and one magnon present. This state has therefore parity (+) since 
Pe = {—),Psw = (+) and Pm = {—)■ We approximate the wave-function as consist- 
ing of two factors existing in different Hilbert spaces. A factor from the end excitations 
and a factor, corresponding to the m-magnon wave-function. From the above dis- 

cussion we see that for a single magnon the magnon part of the wave-function, "^m, be- 
comes \E'i = yj'^jsmkix, rii = 1. From Eq. (H) we then see that < >2+ — < 

>i-= sin^ /cix, which is shown as the solid line in Fig. Here the subtraction 
of < >i- essentially removes any contribution from the end-excitations. In Fig. |I| we 
show < >2+ — < 5'f >!-. An excellent agreement is evident. Also shown in Fig. |I] is 
the local bond energy ef^ =< Sj ■ Sj+i >2+ — < Sj • Sj+i >i-. The dilute boson model 
predicts an energy density, e(x) ^ ^J^JLi^i^ ~ ignoring the Oi^L'"^) kinetic energy. 
Therefore e(x) should be proportional to S^{x) in Eq. (|^, the proportionality factor being 
the gap, A. In Fig. [1| this prediction is shown as the dotted line. From Eq. we can now 
extract values for the gap. A, and the velocity, v. The fit of A2i(L) = E2+{L) — Ei-{L) to 
Eq. d) is excellent, and we obtain A2i(L) = 0.4107(1) + 74.7(4)(L - 1)-^ + 0{[L - 

6 



with = 4.55. We see that A = 0.4107(1), v = 2.49(1). The value of A is in excellent 
agreement with what was previously obtained |1]3. The value of v is in good agreement 
with the value v = 2.46 that can be extracted from exact diagonalization |Q, and the value 



V ~ 2.36 obtained from 1/5* expansions [|T5[. It is also in good agreement with the experi- 
mental results on NENP 0, v ~ 2.45. The coefficient, 74.7(4), in front of the (L - 1)^^ 



term, that determines v, differs marginally from what was obtained by White [IH, (67.9), 
due to the use of a different polynomial form. 

The lowest lying 2-magnon state has parity (— ) since Pe = {—),Psw = {+),Pm = 
(+) and total magnetization M=m+1=3. The magnon part of the wave-function is \E'2 = 
-^^[sin kiXi sin ^22^2 ~ sin ^13^2 sin k2Xi]e{xi, X2) with rii = 1, n2 = 2. Note that under parity 
Xi ^ L — Xi — 1, sin kiXi is even, sin k2Xi is odd and e(xi, X2) is odd resulting in Psw = (+)• 
We now obtain < S'f >3- — < S- >i-= -^;3Y{sin^ fcix + sin^ A;2a;}, which is shown as the 
solid line in Fig. |^. The dotted line represents the theoretical prediction for the local bond 
energy which is also shown in Fig. |^. Again we fit A3i(L) = E^--{L) — Ei-{L) to Eq. (^ 
and we obtain A3i(L) = 0.823(1) + 359(5)(L - 1)"^ + 0{[L - 1]-^). The constant term 
should be 2A in good agreement with the value of A obtained above. Since in this case 
we expect ni = l,n2 = 2 and therefore Z)^? = 5, the coefficient in front of the {L — 1)~^ 
term should be 5 times greater than what we found for the 2~ level. Clearly this is the case: 
we obtain Y^nf = 4.80(6), and the only values of rij consistent with our results are indeed 
rii = l,n2 = 2. 

The three magnon state with M = m + 1 = 4 has parity (+) by the same arguments 
as above. The wave-function now has 6 terms and we obtain < >4+ — < >i- = 
j^{sm^ kix + sin^ k2X + siii?k^x}, with rii = 1, -0.2 = 2, = 3. This expression is shown as 
the solid line in Fig. ^, along with the numerical results for < >4+ — < >i- and the 
local bond energy < Sj ■ Sj+i >4+ — < Sj ■ Sj+i >i-. The dotted line is the prediction for 
the local bond energy. Again good agreement is evident between theory and the numerical 
results is seen. Fitting the energy gap to Eq. (|^) we find that the (L — 1)^^ term now has 
a coefficient of 1030(150). Thus in this case, if we use the value of v determined above. 



^nf = 14(2). This is only consistent with the values rii = l,n2 = 2,723 = 3. The DMRG 
method works progressively worse for states with higher magnetization and the agreement 
between the theory and the numerical results is therefore not as spectacular for this state 
as for the states previously discussed. 

In summary we find that the numerical results are in excellent agreement with Eq. 
with V = 2.49(1), A = 0.4107(1). The magnons behave as bosons with repulsive interactions 
among themselves and with the end excitations. The lowest energy state for given M has 
the Bloch form of Eq. (^, implying the validity of Eq. (^). 
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TABLES 

TABLE I. The spectrum of the L = 100 open = 1 antiferromagnetic Heisenberg chain. 



— E 

1- 138.940086 
2+ 138.522461 
3- 138.08557 
4+ 137.603 
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FIGURES 

FIG. 1. The open circles represent < Sf >2+ — < >i-- The sohd line is the expression 

given in the text. Also shown, by triangles, is < Sj • Sj+i >2+ — < Sj • Sj+i >i-. The dotted line 
is the prediction for this local bond energy. 

FIG. 2. The open circles represent < >3- — < >i-- The solid line is the expression 
given in the text. Also shown, by triangles, is < • Sj+i >3- — < Sj • S^+i >i-. The dotted line 
is the prediction for this local bond energy. 

FIG. 3. The open circles represent < Sf >4+ — < Sf >i-. The solid line is the expression 
given in the text. Also shown, by triangles, is < Sj • S^+i >4+ — < • Sj+i >i-. The dotted line 
is the prediction for this local bond energy. 
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